# Define more focused parameter bounds based on marketing intuition
# We'll make these significantly tighter than before but not directly using the known values
hill_bounds = {
'TV': {
'alpha': (20000, 30000), # Higher impact medium
'gamma': (0.6, 0.8), # Moderate diminishing returns
'k': (30000, 40000) # Higher saturation point
},
'Digital': {
'alpha': (15000, 25000), # Moderate-high impact
'gamma': (0.5, 0.7), # Moderate-high diminishing returns
'k': (20000, 30000) # Moderate saturation point
},
'Radio': {
'alpha': (10000, 15000), # Lower impact
'gamma': (0.7, 0.9), # Lower diminishing returns
'k': (10000, 15000) # Lower saturation point
},
'Print': {
'alpha': (7000, 10000), # Lowest impact
'gamma': (0.8, 1.0), # Lowest diminishing returns
'k': (8000, 12000) # Lowest saturation point
}
}
print("Optimizing Hill function parameters...")
# L2 regularization penalty
l2_penalty = 10
# Optimize Hill parameters for each channel
optimal_hill_params = {}
channel_impacts = {}
# Create time series cross-validation setup - reusing from previous section
from sklearn.model_selection import TimeSeriesSplit
n_splits = 3
tscv = TimeSeriesSplit(n_splits=n_splits)
for channel in channels:
print(f" Processing {channel}...")
# Get adstocked values - assuming they were calculated in previous cell
adstocked_values = data[f"{channel}_Adstocked"].values
# Initialize a randomized search approach to avoid local minima
best_r2 = -float('inf')
best_params = None
num_random_starts = 8 # Try more random starting points for better results
# Define objective function for Hill parameter optimization
def hill_objective(params, train_indices, val_indices=None):
alpha, gamma, k = params
# Apply hill transformation
impact = hill_transform(adstocked_values, alpha, gamma, k)
# If validation indices provided, use them, otherwise use training indices
eval_indices = val_indices if val_indices is not None else train_indices
# Calculate R²
correlation = np.corrcoef(impact[eval_indices],
data['Total_Sales'].values[eval_indices])[0, 1]
r_squared = correlation ** 2
# Add L2 regularization with less emphasis than before
l2_reg_term = (
l2_penalty * 0.01 * ((alpha - 20000)/10000)**2 + # Reduced regularization weight
l2_penalty * 0.01 * ((gamma - 0.7)/0.3)**2 + # Reduced regularization weight
l2_penalty * 0.01 * ((k - 20000)/10000)**2 # Reduced regularization weight
)
# Return negative R² plus regularization (to minimize)
return -r_squared + l2_reg_term
# Try multiple optimization runs with different starting points
for _ in range(num_random_starts):
cv_hill_params = []
cv_hill_scores = []
for train_idx, val_idx in tscv.split(data['Total_Sales']):
# Set bounds for Hill parameters
bounds = [
hill_bounds[channel]['alpha'],
hill_bounds[channel]['gamma'],
hill_bounds[channel]['k']
]
# Start with randomized point within the bounds
alpha_start = np.random.uniform(bounds[0][0], bounds[0][1])
gamma_start = np.random.uniform(bounds[1][0], bounds[1][1])
k_start = np.random.uniform(bounds[2][0], bounds[2][1])
initial_guess = [alpha_start, gamma_start, k_start]
# Define the objective function for the optimizer
obj_func = lambda p: hill_objective(p, train_idx, val_idx)
# Run optimization
from scipy.optimize import minimize
result = minimize(
obj_func,
initial_guess,
bounds=bounds,
method='L-BFGS-B',
options={'maxiter': 200} # More iterations for better convergence
)
cv_hill_params.append(result.x)
cv_hill_scores.append(-result.fun) # Convert back to positive R²
# Calculate average score for this random start
avg_score = np.mean(cv_hill_scores)
if avg_score > best_r2:
best_r2 = avg_score
# Use median of parameters from all CV folds
best_params = [
np.median([p[0] for p in cv_hill_params]),
np.median([p[1] for p in cv_hill_params]),
np.median([p[2] for p in cv_hill_params])
]
# Get the best parameters
optimal_alpha, optimal_gamma, optimal_k = best_params
print(f" Optimal Hill parameters for {channel}:")
print(f" alpha: {optimal_alpha:.2f}")
print(f" gamma: {optimal_gamma:.4f}")
print(f" k: {optimal_k:.2f}")
# Store optimized parameters
optimal_hill_params[channel] = {
'alpha': optimal_alpha,
'gamma': optimal_gamma,
'k': optimal_k
}
# Calculate final impact with optimal parameters
impact = hill_transform(adstocked_values, optimal_alpha, optimal_gamma, optimal_k)
# Store impact and create Hill-transformed variable
channel_impacts[channel] = impact
data[f"{channel}_Hill"] = impact
# Use Ridge regression for Hill-transformed variables
X_hill = data[[f"{channel}_Hill" for channel in channels]]
# Scale the inputs for better numerical stability
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import RidgeCV
from sklearn.metrics import r2_score, mean_squared_error
hill_scaler = StandardScaler()
X_hill_scaled = hill_scaler.fit_transform(X_hill)
# Use cross-validation to find optimal alpha - with more moderate values
alphas = [0.01, 0.1, 1.0, 10, 100, 1000] # Same alphas as in simple ridge model
hill_ridge_cv = RidgeCV(alphas=alphas, cv=5)
hill_ridge_cv.fit(X_hill_scaled, data['Total_Sales'])
# Get optimal alpha
hill_best_alpha = hill_ridge_cv.alpha_
print(f"\nAdstock + Hill Media Model:")
print(f"Optimal regularization strength (alpha): {hill_best_alpha}")
# Get coefficients in original scale
hill_coefs = hill_ridge_cv.coef_ / hill_scaler.scale_
hill_intercept = hill_ridge_cv.intercept_ - np.sum(hill_ridge_cv.coef_ * hill_scaler.mean_ / hill_scaler.scale_)
# Store results in a format similar to the simple ridge model
hill_coef_df = pd.DataFrame({
'variable': ['const'] + [f"{channel}_Hill" for channel in channels],
'coef': [hill_intercept] + list(hill_coefs)
})
print("\nHill Ridge Regression Coefficients:")
print(hill_coef_df)
# Make predictions
data['predicted_hill'] = hill_ridge_cv.predict(X_hill_scaled)
# Calculate model metrics
hill_r2 = r2_score(data['Total_Sales'], data['predicted_hill'])
hill_rmse = np.sqrt(mean_squared_error(data['Total_Sales'], data['predicted_hill']))
print(f"\nAdstock + Hill model R²: {hill_r2:.4f}")
print(f"Adstock + Hill model RMSE: ${hill_rmse:.2f}")
print(f"Improvement over adstock model: {(hill_r2 - adstock_r2) * 100:.2f} percentage points")
print(f"Total improvement over simple model: {(hill_r2 - ridge_r2) * 100:.2f} percentage points")
# Plot comparison of all models - assuming the previous models' predictions exist
plt.figure(figsize=(15, 6))
plt.plot(data['Time_Period'], data['Total_Sales'], marker='o', label='Actual Sales')
plt.plot(data['Time_Period'], data['predicted_simple_ridge'],
linestyle='--', marker='x', alpha=0.5, label='Simple Media Model (Ridge)')
plt.plot(data['Time_Period'], data['predicted_adstock_ridge'],
linestyle='--', marker='+', alpha=0.7, label='Adstock Media Model (Ridge)')
plt.plot(data['Time_Period'], data['predicted_hill'],
linestyle='--', marker='s', label='Adstock + Hill Media Model (Ridge)')
plt.title('Comparison of Media Models: Simple vs Adstock vs Hill (all with Ridge Regression)', fontsize=14)
plt.legend()
plt.xticks(rotation=90)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Added: Plot diminishing returns curves for each channel
plt.figure(figsize=(15, 10))
for i, channel in enumerate(channels):
plt.subplot(2, 2, i+1)
params = optimal_hill_params[channel]
alpha, gamma, k = params['alpha'], params['gamma'], params['k']
# Generate response curve
max_spend = data[f"{channel}_Adstocked"].max()
spend_levels = np.linspace(0, max_spend * 1.5, 100)
response = hill_transform(spend_levels, alpha, gamma, k)
# Plot the curve
plt.plot(spend_levels, response)
# Mark the current max spend point
current_response = hill_transform(max_spend, alpha, gamma, k)
plt.scatter([max_spend], [current_response], color='red', zorder=5)
plt.axvline(x=max_spend, color='gray', linestyle='--', alpha=0.5)
plt.axhline(y=current_response, color='gray', linestyle='--', alpha=0.5)
# Mark 50% response point
plt.axhline(y=alpha/2, color='green', linestyle='--', alpha=0.5)
plt.axvline(x=k, color='green', linestyle='--', alpha=0.5)
plt.title(f'{channel} Diminishing Returns Curve')
plt.xlabel('Adstocked Spend Level')
plt.ylabel('Sales Impact ($)')
plt.grid(True, alpha=0.3)
# Annotate key points
plt.annotate(f'Current Max: ${max_spend:.0f}',
xy=(max_spend, current_response),
xytext=(max_spend*0.7, current_response*0.9),
arrowprops=dict(arrowstyle='->'))
plt.annotate(f'Half-max at k=${k:.0f}',
xy=(k, alpha/2),
xytext=(k*1.2, alpha/2*0.8),
arrowprops=dict(arrowstyle='->'))
plt.tight_layout()
plt.show()
# Display channel diminishing returns analysis
print("\n==== Channel Diminishing Returns Analysis ====")
for channel in channels:
params = optimal_hill_params[channel]
alpha, gamma, k = params['alpha'], params['gamma'], params['k']
print(f"\n{channel} Channel Diminishing Returns:")
print(f" Alpha: {alpha:.2f} - Maximum potential impact")
print(f" Gamma: {gamma:.4f} - Shape parameter (lower = stronger diminishing returns)")
print(f" K: {k:.2f} - Half-saturation constant")
# Generate response curve for visualization
max_spend = data[f"{channel}_Adstocked"].max()
spend_levels = np.linspace(0, max_spend * 1.5, 100)
response = hill_transform(spend_levels, alpha, gamma, k)
# Calculate elasticity at key spend levels (25%, 50%, 75%, 100% of max observed spend)
spend_points = [0.25 * max_spend, 0.5 * max_spend, 0.75 * max_spend, max_spend]
print(f" Marginal Returns Analysis:")
for i, spend in enumerate(spend_points):
# Calculate elasticity (approximate using small delta)
delta = 0.01 * spend
response1 = hill_transform(spend, alpha, gamma, k)
response2 = hill_transform(spend + delta, alpha, gamma, k)
elasticity = ((response2 - response1) / response1) / (delta / spend)
# Calculate percent of maximum impact
percent_of_max = (response1 / alpha) * 100
print(f" At {i+1}0% of max spend (${spend:.2f}):")
print(f" Impact: ${response1:.2f} ({percent_of_max:.1f}% of maximum potential)")
print(f" Elasticity: {elasticity:.4f}")
# Interpretation of diminishing returns
if elasticity > 0.8:
print(f" Interpretation: Strong returns, minimal diminishing effects")
elif elasticity > 0.5:
print(f" Interpretation: Good returns, moderate diminishing effects")
elif elasticity > 0.2:
print(f" Interpretation: Fair returns, significant diminishing effects")
else:
print(f" Interpretation: Poor returns, severe diminishing effects")
